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The  Naval  Space  Command  maintains  a  database  of  orbital  elements  for  the 
objects  in  the  space  catalog.  This  report  is  a  documentation  of  the  software 
which  inputs  new  observations  of  a  satellite  and  its  old  element  set  and  outputs  a 
new  element  set.  Topics  covered  include:  mathematics  of  batch  least  squares 
differential  correction  process,  definition  of  fit  span  and  passes,  calculation  of 
residuals  and  partials,  inclusion  of  historical  data,  solution  to  normal  equations, 
iterations  and  tolerances. 

INTRODUCTION 

The  Naval  Space  Command  maintains  a  database  of  element  sets  for  roughly  9,000 
Earth- orbiting  objects,  which  is  essentially  the  U.S.  Space  Command  satellite  catalog. 
NAVSPACECOM  receives  about  270,000  observations  per  day  and  performs  an  average 
of  18,000  element  set  updates  per  day.  About  98.5%  of  the  element  sets  are  updated, 
without  human  intervention,  by  computer  software  called  AUTODC.  The  purpose  of  our 
report  is  to  elucidate  the  technical  aspects  of  AUTODC  for  the  astrodynamics 
community. 

MATHEMATICAL  BACKGROUND 

AUTODC  uses  a  batch  least  squares  differential  correction  process.  This  is  an 
iterative  method  which  updates  the  orbital  elements  so  that  the  propagated  orbit  fits  the 
observations  in  the  least  squares  sense,  i.e.,  the  sum  of  the  squares  of  the  residuals  is 
minimized.  In  this  section,  we  outline  the  mathematical  steps  in  the  classical  method  for  a 
batch  least  squares  differential  correction  of  an  orbital  element  set;  for  further  details,  see 
Vallado  (Ref.l). 
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(i)  Pick  an  initial  nominal  state 


X 


nominal 


*1 

*8 


These  elements  will  be  updated  with  each  of  the  following  iterations  until  the  differential 
corrections  are  smaller  than  the  prescribed  tolerances. 

(ii)  Compute  the  values  of  the  observed  parameters  Yc  at  N  times  corresponding  to  the 
observations  Y0  (assume  each  observation  set  includes  6  numbers): 
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(iii)  Compute  the  residuals  or  “O-Cs”  (Ya  -  Yc),  (observed  minus  calculated  parameters). 
These  can  be  arranged  to  form  the  6Nxl  column  matrix 
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(iv)  Calculate  the  partial  derivatives 

dYc  =  dYc  8(r,v) 
dX  3(r,  v)  dX 
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at  each  time  corresponding  to  the  observations.  Here  ( n,rj,rK )  and  ( vi,vj,vk )  denote  the 
Cartesian  coordinates  of  the  position  and  velocity  vectors.  These  can  be  arranged  to  form 
the  6Nx8  matrix 
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(v)  Form  the  normal  equations: 

At  W Ax  =  At  W  b 

Here  W  denotes  a  6Nx6N  diagonal  weighting  matrix,  and 

*1 

x  =  : 

*8 

T  T 

are  the  differential  corrections.  Note  that  A  IT  A  is  an  8x8  matrix,  and  that  A  W  b  is  an 
8x1  matrix. 

(vi)  Solve  the  normal  equations: 

x  =  (ArWAJl  ArWb 

(vii)  Update  the  elements: 

^fnew  =  ^flast  +  X 

(For  the  first  iteration,  Flast  is  Xnomina|.) 


(viii)  Apply  tests  to  determine  if  iterations  should  continue.  If  so,  return  to  step  (ii). 

PROGRAMMED  PROCEDURE 


Now  we  outline  aspects  of  the  differential  correction  process  as  implemented  in  the 
NAVSPACECOM  software.  The  software  has  evolved  from  many  years  of  experience 
with  the  maintenance  of  satellite  orbital  elements  where  the  algorithms  must  be  able  to 
also  handle  atypical  cases  which  arise  in  an  operational  environment.  These  can  include 
cases  where  there  is  insufficient  or  poorly  distributed  observation  data  over  the  orbit  of 
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the  satellite,  or  cases  where  observations  have  been  mistagged  to  the  incorrect  satellite. 
The  quality  of  the  available  observation  data  also  affects  how  well  the  DC  process  can 
overcome  these  limitations  in  the  data.  Some  of  the  formulas  and  logic  used  in  the 
software  were  empirically  derived  based  on  experience  with  the  satellite  element  set 
correction  process. 

Documentation  of  these  methods  was  initiated  by  Perini  (Ref.  2). 

Units 

The  AUTODC  software  internally  uses  canonical  units.  The  Earth  radius  R ©  =  ER 
and  gravitational  parameter  ju®  are  taken  to  be  unity.  One  canonical  time  unit  CTU  is  then 
the  time  it  takes  a  hypothetical  satellite  to  travel  one  radian  on  its  way  around  a  circular 
orbit  of  radius  R ©. 

To  relate  these  units  to  metric  values,  the  Earth  radius  and  gravitational  parameter  are 
taken  to  be: 


R®  =  6378. 135  km,  ju®  =  398597.62579588  km3/sec2 
The  value  for  the  canonical  time  unit  is  then  derived  from  the  formula 


CTU  =J—  =  806.813594299989177  sec 

V© 

Time  is  measured  in  canonical  time  units  from  12:00  UTC  on  January  1,  2000.  The 
NAVSPACECOM  inertial  coordinate  system  is  also  associated  with  this  date. 

Leap  seconds  are  not  accounted  for  in  the  conversion  between  external  (wall  clock) 
time  and  internal  (canonical)  time  units.  If  a  leap  second  occurs  within  the  fit  span,  the 
software  adds/subtracts  the  time  tick  to  the  appropriate  observation  times. 

The  software  uses  radians  for  angle  values. 

Observation  Types 

The  collection  of  tracking  assets  available  to  U.S.  Space  Command  for  tracking 
satellites  non-cooperatively  is  known  as  the  Space  Surveillance  Network  (SSN).  This 
"system  of  systems"  presently  consists  of  several  dozen  sensors,  including  mechanical 
radars,  phased  array  radars,  optical  telescopes,  passive  radio  direction-finding  equipment, 
and  the  Naval  Space  Surveillance  System.  The  latter  system  is  commonly  known  as  the 
Fence,  is  operated  by  Naval  Space  Command,  and  supplies  unique  data  types  discussed 
below.  At  present,  AUTODC  can  handle  the  following  combinations  of  data  types  from 
the  SSN: 
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Azimuth  and  elevation  angles 
Right  ascension  and  declination  angles 
Azimuth,  elevation,  and  range 

Earth-fixed  Greenwich  (EFG)  Cartesian  position  (X,Y,Z) 

Azimuth,  elevation,  range,  plus  EFG  sensor  location 

The  first  three  data  types  come  from  sensors  whose  positions  are  in  the  database.  The  last 
data  type  is  specifically  intended  for  mobile  sensors,  such  as  satellite-based  observing 
platforms.  Some  sensors  report  right  ascension  and  declination,  but  these  angles  are 
transformed  into  apparent  azimuth  and  elevation  for  AUTODC  processing.  Some  sensors 
can  also  report  rates  of  change  of  azimuth,  elevation,  or  range,  but  currently  rate  data  are 
not  used  for  updating  cataloged  orbits. 

The  NAVSPACECOM  Fence  provides,  via  specialized  real-time  interferometric 
processing,  a  pair  of  direction  cosines  of  the  apparent  line  of  sight  of  a  satellite  as  it 
passes  through  the  Fence  beam.  The  cosines  are  reckoned  along  nominal  "East-West"  and 
"North-South"  axes  at  each  of  six  receiver  stations  located  along  a  great-circle  arc  across 
the  southern  United  States.  These  local  observation  axes  differ  from  the  true  geographical 
directions  by  a  single  rotation  in  azimuth,  different  for  each  station.  By  triangulating  the 
apparent  lines  of  sight  from  several  stations,  it  is  possible  to  produce  an  EFG  Cartesian 
position  for  each  satellite  pass.  However,  this  triangulated  (X,Y,Z)  data  type  is  prone  to 
systematic  errors,  and  NAVSPACECOM  uses  only  the  single-station  direction  cosines  for 
updating  the  satellite  catalog.  The  real-time  reduction  of  raw  Fence  data  also  provides 
estimates  of  the  direction  cosine  rates,  bistatic  Doppler  shift,  and  Doppler  rate.  The 
Doppler  shift  is  used  to  help  associate  observations  with  known  orbits. 

Element  Database 


Several  different  versions  of  each  element  set  are  kept  on  the  NAVSPACECOM 
database,  for  comparison  purposes: 


initial 

operational 

SSC 


first  element  set  cataloged  at  NAVSPACECOM 

most  recently  updated  NAVSPACECOM  element  set  for  the  object 

most  recent  element  set  received  from  U.S.  Space  Command  for  the 
object  (SSC  was  acronym  for  Space  Surveillance  Center) 


field  =  most  recent  element  set  provided  to  the  Space  Surveillance  Network 
sensors 


The  operational  element  set  update  is  initiated  based  on  the  arrival  of  new  observation 
data  into  the  system.  If  the  new  observation  residuals  calculated  with  the  current 
operational  element  set  exceed  pre-established  tolerances,  an  update  of  the  element  set  is 
initiated.  Or,  if  it  has  been  more  than  23  hours  since  the  last  update,  an  update  is 
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initiated.  However,  if  no  data  is  received  for  a  satellite,  an  automatic  update  of  the 
element  set  is  not  initiated. 

The  field  element  set  update  may  only  be  initiated  when  NAVSPACECOM  is 
activated  for  providing  element  sets  to  the  sensors.  Then,  each  field  set  update  is  initiated 
when  the  position  difference  with  the  orbit  determined  by  the  current  operational  set  is 
greater  than  a  tolerance,  usually  12  km. 

An  element  set  for  each  orbiting  object  contains  six  classical  elements  augmented  by 
two  terms  used  to  account  for  atmospheric  drag: 
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classical 


“  M  " 

mean  anomaly 

n 

mean  motion 

n/2 

first  decay  term 

/i/6 

second  decay  term 

e 

eccentricity 

CO 

argument  of  perigee 

Q. 

right  ascension  of  ascending  node 

cos  i 

cosine  of  inclination 

The  second  decay  term  is  normally  zero  and  is  not  corrected.  Epoch  is  the  time  of  the 
most  recent  observation  in  the  last  DC. 

These  elements  are  currently  propagated  by  the  subroutine  PPT3.  The  older 
subroutine  PPT2  and  the  improved  subroutine  PPT3  have  been  documented  by  Solomon 
(Ref.  3)  and  Schumacher  and  Glover  (Refs.  4-5). 

Fit  Span  and  Passes 

The  "fit  span"  is  the  maximum  length  of  time  over  which  the  observations  are  taken 
for  a  full-batch  differential  correction.  If  the  satellite  period  is  greater  than  or  equal  to  600 
minutes,  fit  span  is  based  on  period  P  =  2 7ttn\  otherwise,  fit  span  is  based  on  rate  of 
change  of  period  P  -  (-  2  7r/n2)ri  : 


Period 

Fit  Span 

P  >  800  minutes 

30  days 

600  <  P  <  800  minutes 

15  days 

P  <  600  minutes: 

P  >  -0.0005  minute s/day 

10  days 

-0.001  <  P  <  -0.0005  minutes/day 

7  days 

-0.01  <  P  <  -0.001  minutes/day 

5  days 

P  <  -0.01  minutes/day 

3  days 
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The  AUTODC  software  may  try  to  expand  this  fit  span  (up  to  a  maximum  of  30  days) 
so  as  to  include  more  observations.  Note  also  that  fit  span,  and  many  other  AUTODC 
parameters,  may  alternatively  be  selectively  adjusted  by  satellite  rather  than  be  set 
automatically  by  the  software. 

The  element  set  is  updated  only  if  there  has  been  at  least  one  observation  since  the 
last  DC  attempt  and  there  are  more  than  5  passes  in  the  fit  span.  Here  a  "pass"  is  defined 
as: 

a.  The  set  of  direction  cosine  observations  from  2  or  more  receiver  sites  when  a  satellite 
passes  through  the  fence  beam. 

b.  The  set  of  observations  provided  by  any  other  SSN  sensor  for  one  satellite  track  over 
the  sensor.  In  this  case,  the  first  observation  in  the  fit  span  is  the  first  pass.  A  new 
observation  is  a  pass  only  if  the  time  since  the  last  observation  is  more  than  10  minutes 
or  if  the  new  observation  is  from  a  different  sensor. 

Tolerances 

If  the  absolute  value  of  any  O-C  residual  is  greater  than  a  tolerance,  the  corresponding 
single  observation  parameter  is  excluded  from  the  iteration. 

Currently  the  tolerances  for  the  direction  cosine,  [azimuth,  elevation,  range],  and 
Earth-fixed  [X,  Y,  Z]  O-Cs  are  all  set  to  the  same  number. 

The  tolerances  for  Doppler  and  range-rate  O-Cs  currently  are  set  to  zero,  which 
effectively  excludes  the  corresponding  observations  from  the  differential  correction 
process.  The  code  for  other  types  of  rate  O-Cs  is  not  in  AUTODC. 

The  initial  tolerance  is  taken  to  be 

TOL  =  2  *  INT[max(2<3  -  1,  1)]  nautical  miles 

Here  INT  denotes  real  number  truncation,  and  a  is  the  mean  semimajor  axis  in  Earth 
radii.  Some  example  values  of  the  initial  tolerance  are  as  follows: 


Period 

Semimaior  Axis 

TOL 

TOL 

120  minutes 

1.264  ER 

2  nm 

3.707  km 

225  minutes 

1.921  ER 

4  nm 

7.413  km 

720  minutes 

4.172  ER 

14  nm 

25.95  km 

1440  minutes 

6.623  ER 

24  nm 

44.48  km 

This  value  for  the  tolerance  is  changed  in  subsequent  DCs,  as  explained  later  in  this 
report. 
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Residuals 


Direction  cosine  and  angular  residuals  and  partials  are  multiplied  by  the  computed 
range  R\  this  scaling  allows  the  use  of  a  common  tolerance  with  units  of  length  for 
residuals  of  all  data  types.  In  addition,  scaling  factors  of  \/sin(elev)  and  cos (elev) ,  where 
elev  is  elevation,  are  introduced  in  two  residuals  only;  this  decreases  the  importance  of 
low-level  East-West  cosine  observations  and  amplifies  azimuth  measurements  at  zenith, 
because  O-Cs  greater  than  the  tolerance  are  excluded.. 

The  differences  between  observed  and  calculated  cosines  are  scaled  as 

OC(l)  =  (East -West  cosine  -Re)  *  R/sin(elev) 

OC(2)  =  (North  -  South  cosine  —  R  •  n)  *  R 


where  hats  denote  unit  vectors.  The  E  and  N  unit  vectors  lie  in  a  plane  tangent  to  the 
local  horizontal,  and  are  reckoned  along  and  normal  to  the  great  circle  arc  formed  by  the 
receiver  station  locations. 

The  differences  between  observed  and  calculated  azimuth  and  elevation  are  scaled  as 

OC( 3)  =  [azimuth  -  arctan(R  ■  e/r  ■  n)]*  R  *  cos  (elev) 

OC(A )  =  [elevation  -  arctan(sin(e/ev)/cos(<?/ev))]*  R 

For  [X,  Y,  Z]  observation  data,  the  O-Cs  are  just  the  differences  between  observed 
and  calculated  positions,  but  expressed  in  the  UVW  coordinate  frame,  where  U  is  the 
direction  from  Earth  center  to  satellite,  V  is  the  transverse  direction  in  the  orbital  plane, 
and  W  is  normal  to  the  orbit  plane. 

The  calculated  quantities  are  propagated  to  the  observation  times  by  the  PPT3  theory. 

Partials  of  Calculated  Observations 

The  partials  of  the  observed  parameters  with  respect  to  position  and  velocity  are 
denoted  in  AUTODC  by 


VO  = 


K 

d(r,v) 


R  times  the  partials  of  the  direction  cosines  with  respect  to  position  r  are 


c 

)(r-e) 

/  yv  yv  \  yv 

R*- 

=  E,~ 

R-E/C 

dr. 

J 

v  >  J 

c. 

)(r-n) 

R*- 

V  / 

=  N:~ 

(r-n)r 

dr. 

J 

\  }  „ 

Here  the  subscript  j  =  1,  2,  3  refers  to  one  of  the  Cartesian  components  (may  be  Earth- 
fixed  or  inertial  coordinate  frame). 

R  times  the  partials  of  azimuth  and  elevation  observations  with  respect  to  r  are 


R* ^-[arctan(R  •  e/r  •  n)]=  [(r  ■  N )ej  -  (r  ■  e]n.  ]/cos2 (elev) 


R*— — (elev)  =  [z 
or, 


(r.z) 


)R: 


COS 


(elev) 


where  Z  is  the  unit  vertical  vector,  so  sin  (elev)  =  R  ■  Z  . 
The  partials  of  range  with  respect  to  r  are 

6rt  ' 


For  [X,  Y,  Z]  observations,  the  partials  with  respect  to  r  are 


drj 


(r-0)=0. 


and  similarly  for  V  and  W . 


Partials  of  Position 

The  partials  of  the  osculating  position  and  velocity  with  respect  to  the  mean  elements 
at  epoch  are  denoted  in  AUTODC  by 

PE=  3(r’y) 

nonsingular 

These  partials  are  calculated  with  respect  to  the  “nonsingular”  elements  at  epoch 
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M  +CO  +  D. 


n 


X 


nonsingular 


ft/  2 

ii/6 

ecos  {co  +  Cl) 
esin  (zy  +  fl) 
sin(z‘/2)cosf2 
sin(z/2)sin  Q 


This  choice  of  partials  implies  that  the  differential  corrections  Xnonsinguiar  will  be  in  terms 
of  nonsingular  elements.  The  nonsingular  elements,  which  are  similar  to  the  well-known 
equinoctial  elements,  are  well-defined  for  zero  eccentricity  or  zero  inclination  orbits; 
however,  for  i  =  n  they  are  ambiguous  since  Q.  is  undefined. 

The  PE  partials  are  calculated  from 


dZ 

nonsingular 

Here  BR  denotes  the  6x6  matrix  of  partials  of  the  osculating  (r,  v)  with  respect  to  the  6 
osculating  elements 


PE  =  BR  *  AR  = 


d(r,v). 


dZ  dX 


M  +CO  +  Z2 


a 

sin(z‘/2)cosr2 
si  n  (z/2)  si  n  Q 
ecos(co  + 12) 
esin((3;-t-i2) 


and  AR  denotes  the  6x8  matrix  of  partials  of  the  osculating  elements  Z  with  respect  to  the 
8  epoch  mean  elements  Xnonsinguiar-  The  AR  partials  are  currently  based  on  the  older  PPT2 
orbital  theory. 

Weights  and  Biases 

The  diagonal  entries  of  the  weighting  matrix  W  currently  used  in  AUTODC  are  either 
obtained  from  the  database  or  set  to  a  standard  weight  of  25xl06. 

Biases  are  loaded  from  the  database  and  are  currently  set  to  zero  for  all  sensors. 
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Formation  of  Normal  Equations 


AUTODC  accumulates  the  effect  of  each  set  of  observations  in  the  8x8  matrices 
£neW  =  Fold  +  PET  *VOr*W*VO*PE 
Gnew  =  Gold  +  PET  *VOT*W*OC 

(For  the  first  set  of  observations,  Eaid  and  G0id  are  identically  zero.)  After  iterating  through 
all  observation  sets,  the  end  result  is  the  normal  equations  in  the  form 

E  x  =  G 


Sequential-Batch  Differential  Correction 

The  matrix  A  WA  from  the  previous  DC  is  stored  in  the  database  and  may  be  used  in 
the  current  DC  to  represent  historical  data.  This  provides  an  option  of  performing  a 
sequential-batch  least  squares  differential  correction  with  only  the  new  observations 
received  since  the  last  successful  DC.  The  mathematical  derivation  of  the  basic  equations 
used  in  AUTODC  is  rather  complicated,  and  significantly  different  from  the  form  given 
by  Vallado,  and  so  is  summarized  here. 

The  least- squares  method  can  be  applied  to  any  linear  system 

A  x~b 

where  typically  the  mxn  system  has  many  more  equations  (m)  than  unknowns  (n)  and  is 
inconsistent.  Because  the  range  of  the  matrix  R(A)  (the  column  space)  is  the  orthogonal 
complement  of  the  null  space  of  the  transpose  N(A  ),  then  multiplication  of  the  equation 
by  the  transpose  always  gives  a  consistent  nxn  system  (though  the  solution  may  not  be 
unique,  if  the  variables  are  correlated) 


Ar  A  x  =  Arb 

This  solution  minimizes  the  residual  r  =  b  -  A  x  by  making  it  orthogonal  to  the  column 

T  T 

space  of  A,  so  that  A  r  =  A  (b  -  A  x)  =  0.  Each  of  the  original  m  equations  can  be 
weighted  differently  using  weights  w,  by  including  a  diagonal  mxm  matrix  W  with  entries 
w\  to  give 


Ar  W  A  x  =  Ar  Wb 

Then  the  residual  r  for  this  weighted  system  satisfies  A  W  r  =  A  W  (b  -  A  x)  =  0. 

For  sequential  batch  least-squares,  suppose  that  the  first  batch  of  data  requires  solving 

A\x~b\ 
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resulting  in  a  (weighted)  least- squares  solution  x\.  Then  a  second  batch  of  data  would 
require  solving 


A2x~b2 

To  utilize  both  sets  of  data  would  require  solving  both  equations  simultaneously,  which 
corresponds  to  the  weighted  least-squares  system 

[(At  WA) i  +  (At  WA)2 \  x  =  (A7  Wb)i  +  ( ATWb)2 

or  using  the  first  solution  xi 

[(Ar  WA) i  +  (A7  WA) 2]  v  =  (A7  W A)\  xi  +  (A7  W b)2 

This  form  corresponds  to  that  given  in  Ref.  1.  Of  course,  to  give  the  older  data  less 
weight,  say  by  a  factor/,  one  can  reduce  the  weights  in  W\,  or  equivalently  scale  the 
entries  of  (A7  W A) by  a  factor/2. 

In  differential  correction,  the  form  above  must  be  modified  because  the  systems  are 
solved  for  changes  in  the  elements  rather  than  the  elements  themselves,  and  because  each 
batch  system  is  iterated  to  account  for  nonlinearity,  and  also  because  each  batch  has  a 
different  epoch.  The  system  to  be  solved  for  each  iteration  on  one  batch  of  data  is 

—  AX  -  A  O 

dx 

where  AX  is  the  change  in  the  elements  and  AO  is  the  observed  quantities  minus  the 
corresponding  quantities  calculated  from  the  current  elements  (O-C’s),  and  the  coefficient 
matrix  dO/dX  is  the  partials  of  the  observed  quantities  with  respect  to  the  elements.  Then 
after  k  iterations,  the  weighted  least-squares  solution  for  the  first  batch  satisfies 

(atWa)^'(x^  -  Xf~1)=(ATwl~1bf-1 

where  the  subscripts  indicate  the  batch,  superscripts  indicate  iteration  number,  Xf  is  the 
elements  found  from  the  first  batch  after  k  iterations,  A*-1  =  r)0,  /r)X  f  - 1  , 
b'r'  =  A  Of-1  =Ox-  Cf”1 ,  with  A,a_1  and  Cf-1  being  calculated  from  the  previous 
iteration’s  elements  Xf-1. 

Then  for  the  next  batch,  the  epoch  shifts  from  l\  to  t2.  The  new  system  after  n 
iterations  is 
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where  the  subscript  indicates  elements  at  the  new  epoch.  To  include  the  old  data  as 
well,  the  old  system  at  the  old  epoch  would  be 


df~'(x"  -Xf~' 


Note  that  these  two  systems  apply  to  the  new  elements  X2  at  different  epochs.  Small 
changes  in  elements  can  be  shifted  approximately  to  the  old  epoch  using  the  linearization 

(AX),  =  -^l(AX),  . 

V  hi  ^  V  J,2 

Then  the  old  system  can  be  written,  using  <f>12  for  dXt)  /dXti  , 


Moving  the  known  quantities  to  the  right-hand  side  and  applying  weighted  least  squares 
gives 


<S>\2(ATWAl~''<S>l2(x’i 


=  d> 


T 
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which  simplifies,  since  the  old  elements  Xf  satisfy  the  old  weighted  least  squares 
equation,  to  give 


Hence,  combining  the  old  data  with  the  new  gives  the  combined  least  squares  system 


(atWa)P  +<tf,(ArW'4"'<l>l2]( 
=  (ATWbt'  +K.(ATWAt'(x;  - 


In  terms  of  the  array  names  used  in  AUTODC,  E  is  (atWa)"2  ' ,  G  is  (ArWb)n~\p  is 
(atWa)1~1  ,  EP  is  d>12 ,  and  GP  is  (xf  -  Xf1 )  ,  so  the  system  becomes 
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(E  +  EPt PEP)  AX  =  (G  +  EPt P  GP ). 


This  explains  the  basic  equations  used  by  AUTODC  to  update  E  and  G. 

Note  for  comparison  purposes  that  this  P  is  the  inverse  of  Vallado’s  covariance 
matrix  P. 

In  the  current  version  of  the  AUTODC  code,  the  option  to  perform  a  sequential-batch 
differential  correction  has  been  bypassed.  The  matrix  A  WA  from  the  previous  DC  is  set 
to  zero,  and  all  observations  in  the  fit  span  are  included  in  the  DC  process.  The 
implementation  of  the  historical  data  option  is  being  reexamined. 

Corrections  in  Classical  Elements 


If  only  one  in  a  pair  of  classical  elements,  either  ( e ,  cd)  or  (O,  cos  i),  should  be 
corrected,  the  AUTODC  code  finds  differential  corrections  Xciassicai  in  terms  of  classical 
elements.  For  this  purpose,  the  partials  in  the  8x8  Jacobian  matrix  are  computed: 


TX  = 


dx 


nonsingular 


dx 


classical 


Two  of  these  partials  contain  a  factor  in  the  denominator  of  sin  i,  which  is  omitted  in  the 
AUTODC  code,  to  prevent  a  singularity  at  zero  inclination. 

The  chain  rule  is  used  to  obtain  the  classical  partials: 

d(r,v)  _  0(r,v) 


PE 


classical 


jjj  ^  X  nonsingular 


dx 


classical 


dx 


nonsingular 


dx 


PE*TX 


classical 


Thus  the  classical  ^classical  and  Gciassicai  matrices  are: 

^classical  =  TXT  *  E  TX 
Gdassical  =  TXT  *  G 

The  solution  to  the  normal  equations  will  be  the  classical  correction  ^classical-  The  chain 
rule  is  used  to  convert  back  to  the  nonsingular  correction: 

-Fionsingular  —  TX  :  -^classical 


Solution  to  Normal  Equations 

The  normal  equations  Ex  =  G  are  solved  by  Gauss-Jordan  elimination  with  full 
pivoting.  The  algorithm  is  identical  to  one  in  the  book  Numerical  Recipes  in  Fortran  by 
Press,  et  al.  (Ref.  6). 
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The  AUTODC  algorithm  contains  the  unusual  feature  that,  when  an  off-diagonal 
element  of  E  is  too  large  relative  to  the  corresponding  diagonal  elements,  i.e.  E~j  > 

(SING )EuEjj ,  then  that  row  (for  even-number  calls)  or  column  (for  odd-number  calls)  is 
inactivated.  The  singularity  indicator  SING  is  initially  set  to  near  one.  Then,  if  the  active 
matrix  is  singular,  the  threshold  SING  is  lowered  by  10%  (to  inactivate  more 
rows/columns)  and  the  solution  is  tried  again.  Eventually,  either  a  solution  is  found,  or 
the  whole  matrix  gets  made  inactive  (flagged  by  a  return  value  SING  =  0),  or  the 
threshold  gets  too  small  (SING  <  0.01)  which  indicates  that  the  whole  system  is  singular. 

Updating  Elements 

Once  the  normal  equations  have  been  solved  for  the  nonsingular  corrections  .rnonsinguiar, 
they  may  be  added  to  the  old  elements  to  form  the  new  elements.  First,  the  classical 
elements  A]asl  classical  from  the  last  iteration  are  converted  to  the  nonsingular  elements 

Iflast  nonsingular- 

Next,  the  nonsingular  elements  are  updated: 


Anew  nonsingular  —  Alast  nonsingular  "t"  E(  1 )  :  -^nonsingular 


Here  V(l)  is  a  relaxation  factor  used  to  scale  down  the  differential  corrections  if  any  is  too 
large.  If  \xj\  >  T\ (xj)  for  any  j  =  1,  . . .,  8,  the  scale  factor  is 


V(l) 


=  min  j=i  g 


where  T\(xj)  are  the  numbers  given  below: 


Tv  (xj) 

T2(xj) 

xi=M  +  co+  Q. 

5xl0'2 

10'5 

X2  =  n 

5xl0'4 

10~6 

*3  =  d/2 

5xl0'6 

10'8 

ac~ 

II 

K 

5xl0"8 

1010 

x$  =  e  cos  (<y+  Q.) 

5xl0'2 

10"5 

X(,  =  e  sin  (&>+  fl) 

5xl0'2 

10'5 

x-j  =  e  sin  (i/2)  cos  £1 

5xl0"2 

10~5 

x%  =  e  sin  (i/2)  sin  Q. 

5xl0'2 

10‘5 
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Finally,  the  new  nonsingular  elements  Xnew  nonsingular  are  converted  to  the  new  classical 
elements  Anew  classical- 

iterations 

The  following  tests  are  applied  to  determine  if  the  iterations  should  continue: 

(i)  If  it  is  the  first  iteration,  continue  iterations. 

(ii)  If  \.\j\  >  T2(xj)  for  any  j  =  1,. .  .,8  ,  continue  iterations.  Here  T2{xj)  are  the  numbers 
given  above. 

>  .1,  continue  iterations. 

Here 

RES=t  t[WT(j)],[OC(j)]f 

i= 1  7=1 

is  the  weighted  sum  of  the  squares  of  the  new  residuals.  The  parameter  V(2)  is  defined  by 


V(2)=ixjV(l)G(j) 

7=1 

If  all  of  the  tests  are  met  in  no  more  than  ITERI  iterations,  the  DC  is  deemed  to  have 
converged.  The  parameter  ITERI  is  set  to  60  for  the  first  set  of  iterations,  and 
subsequently  may  be  set  to  20,  or  set  to  0  as  a  flag  to  compute  statistics. 

Best  Elements 

AUTODC  employs  a  scheme  to  iteratively  tighten  tolerances  and  repeat  the  element 
correction  process  in  order  to  obtain  the  best  element  set  for  the  given  satellite.  The  goal 
is  to  lower  the  RMS  while  retaining  an  "acceptable"  amount  of  data  in  the  DC.  The  final 
outcome  is  the  "best"  set  of  mean  elements.  In  this  process,  the  percentage  of  observation 
data  and  resultant  RMS  used  in  the  DC  are  examined  for  both  the  entire  fit  span  and  for 
the  most  recent  data  received  since  the  last  DC  session. 


An  initial  DC  is  attempted  with  the  tolerance  previously  described  as  the  initial 
tolerance  TOL,  for  a  maximum  of  ITERI  =  60  iterations.  If  more  than  85%  of  the 
observations  in  the  fit  span,  and  80%  of  the  new  observations,  are  used,  the  tolerance  is 
shrunk  to 


TOL  =  STOL  *  max 


jl,  INT 


STOL 
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The  parameter  STOL  is  half  the  initial  tolerance,  and 


RMS 


1 


N  6 


lX\OC(j)£  6  N 


i= i  j= i 


is  the  unweighted  root  mean  square  of  the  residuals.  Here  the  RMS  is  calculated  with 
either  all  residuals  in  the  fit  span,  or  just  the  new  residuals,  whichever  yields  the 
maximum  RMS.  Another  DC  iteration  cycle  is  then  attempted. 

This  process  is  repeated  until  the  new  value  for  TOL  is  less  than  the  tolerance  goal  of 
half  the  initial  tolerance.  As  the  tolerance  is  shrunk,  the  DC  is  considered  acceptable  if 
more  than  85%  of  the  observations  for  the  entire  fit  span,  and  50%  of  the  new 
observations,  are  used  in  the  DC. 

If  only  75-85%  of  the  observations  in  the  fit  span,  and  50%  of  the  new  observations, 
were  used  in  the  last  iteration  cycle,  no  more  attempts  are  made  to  shrink  the  tolerance. 
The  DC  is  terminated  and  the  results  accepted. 

During  this  cycle  of  reducing  the  tolerance  and  attempting  another  DC,  the  elements 
from  each  cycle  are  saved  as  the  best  element  set  if  an  acceptable  amount  of  data  was 
retained  in  the  DC.  If  this  process  of  reducing  the  tolerance  does  not  converge  to  the 
tolerance  goal  with  an  acceptable  amount  of  data,  the  process  is  terminated  and  the  best 
element  set  is  accepted  as  the  new  element  set. 

Provisions  are  made  in  AUTODC  if  the  initial  tolerance  is  too  low  for  an  acceptable 
DC  on  the  first  attempt.  In  this  case,  the  tolerance  is  changed  to  TOL  =  1  ER  and  a  DC  is 
performed  with  zero  iterations  to  compute  statistical  information  only.  Then  a  third  DC  is 
performed  with  ITERI  =  20  and 

max  [STOL,RMS]  | 

STOL 


TOL  =  STOL  *  INT 


The  process  then  continues,  as  described  above,  to  find  the  best  element  set. 

CONCLUSION 

This  has  been  a  brief  summary  of  some  of  the  salient  features  of  the 
NAVSPACECOM  differential  correction  process.  Lor  more  explanation,  and  for  the 
Lortran  source  code,  the  reader  is  referred  to  the  documentation  of  Danielson  and 
Canright  (Ref.  7). 

Now  that  this  documentation  has  been  completed  the  next  step  is  to  study  where 
improvements  can  be  made  in  the  software.  Progress  in  this  direction  has  been  made  by 
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Marshall  (Ref.  8),  who  showed  the  benefits  of  Singular  Value  Decomposition  when 
solving  the  normal  equations.  Goals  of  the  improvement  process  include 

(i)  Enabling  the  software  to  update  100%  of  the  element  sets  automatically 

(ii)  Improving  the  accuracy  of  the  state  and  covariance 

(iii)  Decreasing  computer  run  time. 
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